# /// script
# requires-python = ">=3.10"
# dependencies = [
#     "numpy",
# ]
# ///

from pathlib import Path

import numpy as np

FILE = Path("CQ_Tau_70pix_87modes_100samples_moments_new.npz")
FILE_CONV = Path("CQ_Tau_70pix_87modes_100samples_conv_moments_new.npz")
# FILE = Path("MWC758_120pix_150modes_100samples_moments_new.npz")
# FILE_CONV = Path("MWC758_120pix_150modes_100samples_conv_moments_new.npz")
# FILE = Path("HD_135344B_80pix_100modes_100samples_moments_new.npz")
# FILE_CONV = Path("HD_135344B_80pix_100modes_100samples_conv_moments_new.npz")


def main():
    # Read moment samples files
    moments_samples = np.load(FILE)
    moments_samples_conv = np.load(FILE_CONV)
    # Get peak, velocity, line widths
    I = moments_samples["peak"]
    v = moments_samples["centroid"]
    w = moments_samples["width"]
    I_conv = moments_samples_conv["peak"]
    v_conv = moments_samples_conv["centroid"]
    w_conv = moments_samples_conv["width"]

    # Posterior mean and std devs simply with
    I_mean = np.nanmean(I, axis=0)
    I_std = np.nanstd(I, axis=0)


if __name__ == "__main__":
    main()
